home *** CD-ROM | disk | FTP | other *** search
- /*
- File: meshsurf.c
- Authors: Charles Loop
- Last Modified: June 17 1989
- Purpose: Better implementation (faster, uses less memory)
- */
- #include <stdio.h>
- #include "Math.h"
- #include "4D.h"
- #include "mesh.h"
-
- extern Plane *PlaneEquation();
-
-
- /********** DECasteljau stack variables and constants **********/
-
- #define MAX_DEGREE 6
- #define STACK_SIZE 84 /* choose(MAX_DEGREE + 3, 3) = (MAX_DEGREE+1)*(MAX_DEGREE+2)*(MAX_DEGREE+3)/6 */
-
- static PointType4D B[STACK_SIZE];
-
- static int degree;
- static int top;
-
-
- int LOC(h,i,j,k) /* index function */
- int h,i,j,k;
- {
- return top - (k + (h-i)*(h-i+1)/2 + h*(h+1)*(h+2)/6);
- }
-
- /***************************************************************/
-
-
- /************ surface tessellation *****************************/
-
- static Vertex *Q;
-
- static int EDGE_DIVISIONS;
- static int NUMBER_SAMPLES;
-
- int ord(i,j,k) /* index function */
- int i,j,k;
- {
- return NUMBER_SAMPLES - (EDGE_DIVISIONS-i)*(EDGE_DIVISIONS-i+1)/2 - k - 1;
- }
-
- /***************************************************************/
-
-
- DECasteljau(r,s,t)
- double r,s,t;
- {
-
- int h, i, j, k;
- int l0, l1, l2, l3;
-
- for (h = degree-1; h >= 0; h--) {
- for (i = 0; i <= h; i++) {
- for (j = 0; j <= (h-i); j++) {
- k = h - i - j;
-
- l0 = LOC(h,i,j,k);
- l1 = LOC(h+1,i+1,j,k);
- l2 = LOC(h+1,i,j+1,k);
- l3 = LOC(h+1,i,j,k+1);
-
- B[l0].P[X] = r*B[l1].P[X] + s*B[l2].P[X] + t*B[l3].P[X];
- B[l0].P[Y] = r*B[l1].P[Y] + s*B[l2].P[Y] + t*B[l3].P[Y];
- B[l0].P[Z] = r*B[l1].P[Z] + s*B[l2].P[Z] + t*B[l3].P[Z];
- B[l0].P[W] = r*B[l1].P[W] + s*B[l2].P[W] + t*B[l3].P[W];
-
- }
- }
- }
- }
-
- /***************************************************************/
-
- Go(node, r1, r2, r3, s1, s2, s3, t1, t2, t3, l)
- Mesh *node;
- unsigned r1, r2, r3, s1, s2, s3, t1, t2, t3;
- int l;
- {
- Mesh *submesh;
- int i;
-
- if (l) {
- node->type = SUBMESH;
-
- submesh = (Mesh *)calloc(4, sizeof(Mesh));
-
- Go(submesh , (s1+t1)/2, (s2+t2)/2, (s3+t3)/2,
- (t1+r1)/2, (t2+r2)/2, (t3+r3)/2,
- (r1+s1)/2, (r2+s2)/2, (r3+s3)/2, l-1);
-
- Go(submesh+1, r1, r2, r3,
- (r1+s1)/2, (r2+s2)/2, (r3+s3)/2,
- (t1+r1)/2, (t2+r2)/2, (t3+r3)/2, l-1);
-
- Go(submesh+2, (r1+s1)/2, (r2+s2)/2, (r3+s3)/2,
- s1, s2, s3,
- (s1+t1)/2, (s2+t2)/2, (s3+t3)/2, l-1);
-
- Go(submesh+3, (t1+r1)/2, (t2+r2)/2, (t3+r3)/2,
- (s1+t1)/2, (s2+t2)/2, (s3+t3)/2,
- t1, t2, t3, l-1);
-
- CombindBoxes(node->box, submesh, submesh+1, submesh+2, submesh+3);
-
- for (i = 0; i < 3; i++)
- submesh[i].next = submesh+i+1;
-
- node->sub.m = submesh;
-
- }
- else {
- node->type = TRIANGLE;
-
- node->sub.t = (Triangle *)calloc(1,sizeof(Triangle));
-
- node->sub.t->V[0] = &Q[ord(r1, r2, r3)];
- node->sub.t->V[1] = &Q[ord(s1, s2, s3)];
- node->sub.t->V[2] = &Q[ord(t1, t2, t3)];
-
- minmax(node->sub.t->V[0]->x,node->sub.t->V[1]->x,node->sub.t->V[2]->x,
- &node->box[X][MIN],&node->box[X][MAX]);
- minmax(node->sub.t->V[0]->y,node->sub.t->V[1]->y,node->sub.t->V[2]->y,
- &node->box[Y][MIN],&node->box[Y][MAX]);
- minmax(node->sub.t->V[0]->z,node->sub.t->V[1]->z,node->sub.t->V[2]->z,
- &node->box[Z][MIN],&node->box[Z][MAX]);
-
- node->sub.t->pl = PlaneEquation(node->sub.t->V[0],
- node->sub.t->V[1],
- node->sub.t->V[2]);
-
- }
- }
-
- /****************************************************************************/
-
-
- MeshSurf(m, d, level)
- Mesh *m;
- int d, level;
- {
- Mesh *mp;
- Edge *ep;
- int i,j,k,l,I;
- double r, s, tmp;
-
- if (m) {
- if ((m->type == MESH) && (m->sub.m->type == FACE)) {
-
- mp = m->sub.m;
-
- if (mp->sub.e) {
-
-
- if (level) {
-
- degree = d;
- top = (degree+1)*(degree+2)*(degree+3)/6 - 1;
-
- for (i = degree; i > 0; i--)
- for (j = degree-i; j >= 0; j--) {
- k = degree-i-j;
-
- B[LOC(degree,i,j,k)] = *mp->sub.e->p;
- B[LOC(degree,i-1,j+1,k)] = *mp->sub.e->next->p;
- B[LOC(degree,i-1,j,k+1)] = *mp->sub.e->next->next->p;
-
- mp = mp->next;
- }
-
- EDGE_DIVISIONS = 1 << level;
- NUMBER_SAMPLES = (EDGE_DIVISIONS+1)*(EDGE_DIVISIONS+2)/2;
-
- Q = (Vertex *)calloc(NUMBER_SAMPLES,sizeof(Vertex));
-
- I = 0;
- for (i = 0; i <= EDGE_DIVISIONS; i++) {
- r = (double)i/EDGE_DIVISIONS;
- for (j = 0; j <= (EDGE_DIVISIONS-i); j++) {
- s = (double)j/EDGE_DIVISIONS;
-
- DECasteljau(r, s, 1-r-s);
-
- for (l = top; l > top-4; l--)
- for (k = 0; k < 4; k++)
- B[l].P[k] /= B[l].P[W];
-
- Q[I].x = B[top].P[X];
- Q[I].y = B[top].P[Y];
- Q[I].z = B[top].P[Z];
-
- for (k = 0; k < 4; k++) {
- B[top-1].P[k] = B[top-1].P[k] - B[top-3].P[k];
- B[top-2].P[k] = B[top-2].P[k] - B[top-3].P[k];
- }
-
- Q[I].a = B[top-1].P[Y]*B[top-2].P[Z] - B[top-2].P[Y]*B[top-1].P[Z];
- Q[I].b = B[top-1].P[Z]*B[top-2].P[X] - B[top-2].P[Z]*B[top-1].P[X];
- Q[I].c = B[top-1].P[X]*B[top-2].P[Y] - B[top-2].P[X]*B[top-1].P[Y];
-
- tmp = sqrt(Q[I].a*Q[I].a
- + Q[I].b*Q[I].b
- + Q[I].c*Q[I].c);
-
- Q[I].a = Q[I].a/tmp;
- Q[I].b = Q[I].b/tmp;
- Q[I].c = Q[I].c/tmp;
-
- I++;
-
- }
- }
-
- Go(m, EDGE_DIVISIONS, 0, 0, 0, EDGE_DIVISIONS, 0, 0, 0, EDGE_DIVISIONS, level);
-
- }
- }
- } else if (m->type == MESH) {
-
- for (mp = m->sub.m; mp != (Mesh *) 0; mp = mp->next)
- MeshSurf(mp, d, level);
-
- for (i = 0; i < 3; i++) {
- m->box[i][MIN] = m->sub.m->box[i][MIN];
- m->box[i][MAX] = m->sub.m->box[i][MAX];
- }
-
- for (mp = m->sub.m->next; mp != (Mesh *) 0; mp = mp->next)
- for (i = 0; i < 3; i++) {
- if (mp->box[i][MIN] < m->box[i][MIN])
- m->box[i][MIN] = mp->box[i][MIN];
- if (mp->box[i][MAX] > m->box[i][MAX])
- m->box[i][MAX] = mp->box[i][MAX];
- }
-
- }
- }
- }
-
-
-
-